library(rdrobust)
library(rdd)
library(tidyverse)


load("RD_parties_and_violence.RData")

###################################################################
###################################################################
#       Why Programmatic Parties Reduce Criminal Violence:
#                 Theory and Evidence from Brazil
#
#          Camilo Nieto-Matiz                Natan Skigin
#       camilo.nieto-matiz@utsa.edu          nskigin@nd.edu
###################################################################
###################################################################




#================================================
#                   Figure 1 (a)
# Effect of programmatic party on homicide rates
#================================================

# numerical results
summary(rdrobust(bra_all$hom_avg,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd", covs=NULL))

# set wd
prog.rd <- subset(bra_all)
# make the binned averages, first the bins to the left
count <- 1
bin.size <- .02
binx.left <- vector(length=length(seq(-.2, 0, bin.size)))
biny.left <- vector(length=length(binx.left))
last <- -.02
for(j in seq(-2, 0, bin.size)) {
  biny.left[count] <- mean(prog.rd$hom_avg[prog.rd$prog >= j-bin.size & prog.rd$prog < j], na.rm=T)
  binx.left[count] <- (j+last)/2
  last <- j
  count <- count + 1
}

# now the bins on the right
count <- 1
bin.size <- .02
binx.right <- vector(length=length(seq(0,.2, bin.size)))
biny.right <- vector(length=length(binx.right))
last <- 0.2
for(j in seq(0,2, bin.size)) {
  biny.right[count] <- mean(prog.rd$hom_avg[prog.rd$prog >= j-bin.size & prog.rd$prog < j], na.rm=T)
  binx.right[count] <- (j+last)/2
  last <- j
  count <- count + 1
}


temp1 <- data.frame(cbind(binx.left, biny.left))
temp2 <- data.frame(cbind(binx.right, biny.right))
colnames(temp2) <- c("binx.left", "biny.left")
temp <- rbind(temp1, temp2)

prog.rd<- select(prog.rd, prog, hom_avg) %>% dplyr::rename(binx.left=prog, biny.left=hom_avg)
 


ggplot(data=temp, aes(x=binx.left, y =  biny.left)) +
  geom_point(data=subset(temp, binx.left >= 0&binx.left <= .2),color="royalblue2",alpha = .4, size = 4.5) +
  geom_point(data=subset(temp, binx.left <= 0&binx.left >= -.2),color="tomato2",alpha = .4, size = 4.5) +
  geom_vline(xintercept = 0) +
  geom_smooth(data=subset(prog.rd, binx.left >= 0&binx.left <= .2),method = "loess",
              formula = y ~ poly(x, 2), se = TRUE, size = 1.5, fill="blue", color = "gray10", span = 1, alpha = .4)  +
  geom_smooth(data=subset(prog.rd, binx.left < 0&binx.left >= -.2),method = "loess",
              formula = y ~ poly(x, 2),se = TRUE, size = 1.5, fill="red", color = "gray10", span = 1, alpha = .4)  +
  scale_x_continuous("Vote Margin", expand = c(0, 0)) +
  scale_y_continuous("Homicide Rates", expand = c(0, 0)) +
  coord_cartesian(ylim = c(5,25), xlim = c(-.203,.203))+
  theme_bw() + theme(legend.position = "none")  +
  theme(plot.title = element_text(hjust = 0.5, size=25,face="bold"), axis.title.x = element_text(size=24,face="bold")) +
  theme(plot.title = element_text(hjust = 0.5, size=25,face="bold"), axis.title=element_text(size=25,face="bold")) +
  theme(plot.title = element_text(hjust = 0.5, size=25,face="bold"), axis.text.x.bottom =element_text(size=15,face="bold"), axis.text.y.left =element_text(size=15,face="bold")) +
  theme(plot.title = element_text(hjust = 0.5, size=40,face="bold"), title =element_text(size=5,face="bold"))
 
ggsave(file = "paper_figures/Figure_1a.pdf",
       width = 22, height = 12.5, dpi = 600)


#================================================
#                   Figure 1 (b)
# Effect of non-programmatic party on homicide rates
#================================================

# numerical results
summary(rdrobust(bra_all$hom_avg,bra_all$nonprog,c=0,p=1,  
                 all = T, vce = "hc0", kernel="tri",
                 bwselect = "mserd", covs=NULL))

nonprog.rd <- subset(bra_all)
# make the binned averages, first the bins to the left
count <- 1
bin.size <- .02
binx.left <- vector(length=length(seq(-.2, 0, bin.size)))
biny.left <- vector(length=length(binx.left))
last <- -.02
for(j in seq(-2, 0, bin.size)) {
  biny.left[count] <- mean(nonprog.rd$hom_avg[nonprog.rd$nonprog >= j-bin.size & nonprog.rd$nonprog < j], na.rm=T)
  binx.left[count] <- (j+last)/2
  last <- j
  count <- count + 1
}

# now the bins on the right
count <- 1
bin.size <- .02
binx.right <- vector(length=length(seq(0,.2, bin.size)))
biny.right <- vector(length=length(binx.right))
last <- 0.2
for(j in seq(0,2, bin.size)) {
  biny.right[count] <- mean(nonprog.rd$hom_avg[nonprog.rd$nonprog >= j-bin.size & nonprog.rd$nonprog < j], na.rm=T)
  binx.right[count] <- (j+last)/2
  last <- j
  count <- count + 1
}


temp1 <- data.frame(cbind(binx.left, biny.left))
temp2 <- data.frame(cbind(binx.right, biny.right))
colnames(temp2) <- c("binx.left", "biny.left")
temp <- rbind(temp1, temp2)

nonprog.rd<- select(nonprog.rd, nonprog, hom_avg) %>% dplyr::rename(binx.left=nonprog, biny.left=hom_avg)


ggplot(data=temp, aes(x=binx.left, y =  biny.left)) +
  geom_point(data=subset(temp, binx.left >= 0&binx.left <= .2),color="royalblue2",alpha = .4, size = 4.5) +
  geom_point(data=subset(temp, binx.left <= 0&binx.left >= -.2),color="tomato2",alpha = .4, size = 4.5) +
   geom_vline(xintercept = 0) +
  geom_smooth(data=subset(nonprog.rd, binx.left > 0&binx.left <= .2),method = "loess",formula = y ~ poly(x, 2), se = TRUE, size = 1.5, fill="blue", color = "gray10", span = 1, alpha = .4)  +
  geom_smooth(data=subset(nonprog.rd, binx.left <= 0&binx.left >= -.2),method = "loess",formula = y ~ poly(x, 2),se = TRUE, size = 1.5, fill="red", color = "gray10", span = 1, alpha = .4)  +
  scale_x_continuous("Vote Margin", expand = c(0, 0)) +
  scale_y_continuous("Homicide Rates", expand = c(0, 0)) +
  coord_cartesian(ylim = c(5,25), xlim = c(-.203,.203))+
  theme_bw() + theme(legend.position = "none")  +
  theme(plot.title = element_text(hjust = 0.5, size=25,face="bold"), axis.title.x = element_text(size=24,face="bold")) +
  theme(plot.title = element_text(hjust = 0.5, size=25,face="bold"), axis.title=element_text(size=25,face="bold")) +
  theme(plot.title = element_text(hjust = 0.5, size=25,face="bold"), axis.text.x.bottom =element_text(size=15,face="bold"), axis.text.y.left =element_text(size=15,face="bold")) +
  theme(plot.title = element_text(hjust = 0.5, size=40,face="bold"), title =element_text(size=5,face="bold"))

ggsave(file = "paper_figures/Figure_1b.pdf",
       width = 22, height = 12.5, dpi = 600)



#================================================
#                Figure 2 (a)
#  Effect of programmatic party on yearly rates
#================================================

## homicides  - programmatic
summary(m1_1<-rdrobust(bra_all$homlag,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd", covs=NULL))
summary(m2_1<-rdrobust(bra_all$hom0,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd", covs=NULL))
summary(m3_1<-rdrobust(bra_all$hom1,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd", covs=NULL))
summary(m4_1<-rdrobust(bra_all$hom2,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd", covs=NULL))
summary(m5_1<-rdrobust(bra_all$hom3,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd", covs=NULL))
summary(m6_1<-rdrobust(bra_all$hom4,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd", covs=NULL))



# store estimates
spatial.Frame <- data.frame(Coefficient = c(m1_1$coef[3], m2_1$coef[3],m3_1$coef[3], m4_1$coef[3],m5_1$coef[3],m6_1$coef[3]),
                            SE = c(m1_1$se[3], m2_1$se[3],m3_1$se[3],m4_1$se[3],m5_1$se[3], m6_1$se[3]))
spatial.Frame$Year <-rep(c("t-1","t","t+1","t+2", "t+3", "t+4"),1)
spatial.Frame$Var <-  factor(c(rep("PT-PPS-PV-PSDB", 6)))
spatial.Frame$n <- rep(c(1:6),1)

# plot 
spatial.Frame %>%
  mutate(Year=reorder(Year, n)) %>%
  ggplot()+ 
  geom_errorbar(aes(x = Year, ymin = Coefficient - SE*1.96,
                    ymax = Coefficient + SE*1.96),
                width=0, lwd=.7, position = position_dodge(width = .45), color="gray50") +
  geom_errorbar(aes(x = Year, ymin = Coefficient - SE*1.64,
                    ymax = Coefficient + SE*1.64),
                width=0, lwd=1.5, position = position_dodge(width = .45), color="gray10") +
  geom_point(aes(x = Year, y = Coefficient),
             lwd =2,position = position_dodge(width = .45), color="grey10") +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + theme_bw() + ylim(-10,10)+
  theme(legend.position = "bottom",
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        legend.title=element_blank())
 
ggsave(file = "paper_figures/Figure_2a.pdf",
       width = 22, height = 12.5, dpi = 600)


#================================================
#                Figure 2 (b)
#  Effect of non-programmatic party on yearly rates
#================================================

## homicides  - non programmatic
summary(npm1_1<-rdrobust(bra_all$homlag,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd", covs=NULL))
summary(npm2_1<-rdrobust(bra_all$hom0,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd", covs=NULL))
summary(npm3_1<-rdrobust(bra_all$hom1,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd", covs=NULL))
summary(npm4_1<-rdrobust(bra_all$hom2,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd", covs=NULL))
summary(npm5_1<-rdrobust(bra_all$hom3,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd", covs=NULL))
summary(npm6_1<-rdrobust(bra_all$hom4,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd", covs=NULL))



## store estimates
spatial.Framenp <- data.frame(Coefficient = c(npm1_1$coef[3], npm2_1$coef[3],npm3_1$coef[3], npm4_1$coef[3],npm5_1$coef[3],npm6_1$coef[3]),
                              SE = c(npm1_1$se[3], npm2_1$se[3],npm3_1$se[3],npm4_1$se[3],npm5_1$se[3], npm6_1$se[3]))
spatial.Framenp$Year <-rep(c("t-1","t","t+1","t+2", "t+3", "t+4"),1)
spatial.Framenp$Var <-  factor(c(rep("PP-PL-PTB-PSB", 6)))
spatial.Framenp$n <- rep(c(1:6),1)



# plot 
spatial.Framenp %>%
  mutate(Year=reorder(Year, n)) %>%
  ggplot()+ 
  geom_errorbar(aes(x = Year, ymin = Coefficient - SE*1.96,
                    ymax = Coefficient + SE*1.96),
                width=0, lwd=.7, position = position_dodge(width = .45), color="gray50") +
  geom_errorbar(aes(x = Year, ymin = Coefficient - SE*1.64,
                    ymax = Coefficient + SE*1.64),
                width=0, lwd=1.5, position = position_dodge(width = .45), color="gray10") +
  geom_point(aes(x = Year, y = Coefficient),
             lwd =2,position = position_dodge(width = .45), color="grey10") +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + theme_bw() + ylim(-10,20)+
  theme(legend.position = "bottom",
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        legend.title=element_blank())

ggsave(file = "paper_figures/Figure_2b.pdf",
       width = 22, height = 12.5, dpi = 600)

 